#include "AnaDef.h"
#include "XGLUtils.h"

Double_t Gaus2(const Double_t *xx, const Double_t *par)
{
  return XGLUtils::GausN(xx, par, 2, XGLUtils::RawGaus, 3);
}

void demoRepulsive()
{
  TList *ll=new TList;

  const Double_t mean1 = 3;
  const Double_t mean2 = 7;
  TRandom3 ran(123);
  const Int_t nn=10000;
  XGLUtils::FFunc ffunc = Gaus2;
  const Int_t npar = 6;
  TGraph *g1mean=new TGraph; g1mean->SetName("g1mean"); ll->Add(g1mean);
  TGraph *g2mean=new TGraph; g2mean->SetName("g2mean"); ll->Add(g2mean);

  TGraph *g1amp=new TGraph; g1amp->SetName("g1amp"); ll->Add(g1amp);
  TGraph *g2amp=new TGraph; g2amp->SetName("g2amp"); ll->Add(g2amp);

  //sigma 0.1 - 3
  for(Int_t igr=0; igr<30; igr++){
    const Double_t sigma = (igr+1)/10.+0.5;

    const Double_t xmin = -5;
    const Double_t xmax = 15;
    TH1D *hh = new TH1D(Form("hh%d", igr),"",400, xmin, xmax); ll->Add(hh);

    for(Int_t ii=0; ii<nn; ii++){
      hh->Fill(ran.Gaus(mean1, sigma));
      hh->Fill(ran.Gaus(mean2, sigma));
    }

    //hh->Draw();
    Double_t pars[]={hh->GetEntries()/2*hh->GetBinWidth(1), mean1-1, sigma, hh->GetEntries()/2*hh->GetBinWidth(1), mean2+1, sigma};

    const Bool_t kfail = XGLUtils::BinnedLikelihoodFit(hh, ffunc, npar, pars);
    TH1D * hfit = XGLUtils::GetHfit(Form("hfit%d", igr), ffunc, pars, xmin, xmax);ll->Add(hfit);
    TH1D * g1fit = XGLUtils::GetHfit(Form("g1fit%d", igr), XGLUtils::RawGaus, pars, xmin, xmax);ll->Add(g1fit);
    TH1D * g2fit = XGLUtils::GetHfit(Form("g2fit%d", igr), XGLUtils::RawGaus, &(pars[3]), xmin, xmax);ll->Add(g2fit);

    if(!kfail){
      const Int_t np = g1mean->GetN();
      const Double_t xx = 2*sigma/(mean2-mean1);
      g1mean->SetPoint(np, xx, pars[1]/mean1);
      g2mean->SetPoint(np, xx, pars[4]/mean2);

      g1amp->SetPoint(np, xx, pars[0]/(nn*hh->GetBinWidth(1)));
      g2amp->SetPoint(np, xx, pars[3]/(nn*hh->GetBinWidth(1)));

      g1mean->SetMarkerStyle(20);
      g1amp->SetMarkerStyle(26);
      g1mean->SetMarkerColor(kRed);
      g1amp->SetMarkerColor(kRed);
      g1mean->SetLineColor(kRed);
      g1amp->SetLineColor(kRed);

      g2mean->SetMarkerStyle(20);
      g2amp->SetMarkerStyle(26);
      g2mean->SetMarkerColor(kBlue);
      g2amp->SetMarkerColor(kBlue);
      g2mean->SetLineColor(kBlue);
      g2amp->SetLineColor(kBlue);


    }
  }

  TFile *fout=new TFile("outDemoRepul.root","recreate");
  ll->Write();
  fout->Save();
  fout->Close();

  delete fout;
}

int main()
{

  demoRepulsive();

  return 0;
}
